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Abstract— A family of the Block Matching 3-D (BM3D) algo- 
rithms for various imaging problems has been recently proposed 
within the framework of nonlocal patch-wise image modeling 1 1 1, 
|2|. In this paper we construct analysis and synthesis frames, 
formalizing the BM3D image modeling and use these frames 
to develop novel iterative deblurring algorithms. We consider 
two different formulations of the deblurring problem: one given 
by minimization of the single objective function and another 
based on the Nash equilibrium balance of two objective functions. 
The latter results in an algorithm where the denoising and 
deblurring operations are decoupled. The convergence of the 
developed algorithms is proved. Simulation experiments show 
that the decoupled algorithm derived from the Nash equilibrium 
formulation demonstrates the best numerical and visual results 
and shows superiority with respect to the state of the art in 
the field, confirming a valuable potential of BM3D-frames as an 
advanced image modeling tool. 

I. Introduction 

WE consider image restoration from a blurry and noisy 
observation. Assuming a circular shift-invariant blur 
operator and additive zero-mean white Gaussian noise the 
conventional observation model is expressed as 

z = Ay + cre, (1) 

where z, y 6 R N are vectors representing the observed 
and true image, respectively, A is an N x N blur matrix, 
e ~ JV(Onxi, Inxn) is a vector of i.i.d. Gaussian random 
variables, and a is the standard deviation of the noise. The 
deblurring problem is to reconstruct y from the observation 
z. The most popular approach is to formulate reconstruction 
as a variational optimization problem, where the desired so- 
lution minimizes a criterion composed of fidelity and penalty 
terms. The fidelity ensures that the solution agrees with the 
observation, while the penalty provides regularization of the 
optimization problem through a prior image model. Typically, 
the fidelity term is derived from the negative log-likelihood 
function. For the Gaussian observation model (HJ the fidelity 

term has the form — 2 ||z — Ay|| 2 , and the minimization 
criterion is given as 

J=^\\z-Ayf 2 +T-pen{y), (2) 

where ||-|| 2 stands for the Euclidean norm, pen(-) is a penalty 
functional and r > is a regularization parameter. 
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Image modeling lies at the core of image reconstruction 
problems. Recent trends are concentrated on sparse represen- 
tation techniques, where the image is assumed to be defined 
as a combination of few atomic functions taken from a certain 
dictionary. It follows that the image can be parameterized 
and approximated locally or nonlocally by these functions. 
To enable sparse approximations, the dictionary should be 
rich enough to grasp all variety of the images. Clearly, bases 
are too limited for this task and one needs to consider 
overcomplete systems with a number of elements essentially 
larger than the dimensionality of the approximated images. 
Frames are generalization of the concept of basis to the case 
when the atomic functions are linearly dependent and form an 
overcomplete system (3). There is a vast amount of literature 
devoted to the sparsity based models and methods for imaging. 
An excellent introduction and overview of this area can be 
found in the recent book [4 |. 

The contribution of this paper concerns three main aspects 
of image deblurring: image modeling, variational problem 
formulation, and algorithmic reconstruction. 

First, the BM3D image modeling developed in [ 1 1 is formal- 
ized in terms of the overcomplete sparse frame representation. 
We construct analysis and synthesis BM3D-frames and study 
their properties. The analysis and synthesis developed in 
BM3D are interpreted as a general sparse image modeling 
applicable to variational formulations of various image pro- 
cessing problems. 

Second, we consider two different formulations of the 
image deblurring problem: one given by minimization of the 
objective function and another based on the Nash equilibrium. 
The latter approach results in an algorithm where the denoising 
and the deblurring operations are decoupled. 

Third, it is shown by simulation experiments that the best 
image reconstruction both visually and numerically is obtained 
by the algorithm based on decoupling of blur inverse and noise 
filtering. To the best of our knowledge, this algorithm provides 
results which are the state-of-art in the field. 

Here we extend and develop our preliminary ideas sketched 
in 0. The BM3D frames are now constructed explicitly, 
taking into account the particular form of the 3D transform. 
Proofs of the frame properties are presented. We develop algo- 
rithms for the analysis and synthesis-based problem formula- 
tions introduced in [5 1 and provide their convergence analysis. 
The problem formulation based on the Nash equilibrium and 
the corresponding decoupled deblurring algorithm are novel 
developments. 

The paper is organized as follows. We start from a presen- 
tation of the BM3D image modeling and introduce BM3D- 
frames (Section Hill. The variational image reconstruction is a 
subject of Section [Hi] The algorithms based on the analysis 
and synthesis formulations are derived in this section. The al- 
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gorithm based on the Nash equilibrium is presented in Section 
IIVI Convergence results for the proposed algorithms are given 
in Section [V] Implementation of the algorithms is discussed in 
Section[VT] The experiments and comparison of the algorithms 
are given in Section IVHI In Section IVIIII we discuss the 
principal differences of the decoupled formulation compared to 
the analysis and synthesis formulations. Concluding remarks 
are done in the last section. Proofs of mathematical statements 
are given in Appendix. 

II. OVERCOMPLETE BM3D IMAGE MODELING 

BM3D is a nonlocal image modelling technique based on 
adaptive, high order groupwise models. Its detailed discussion 
can be found in [6|. Below, using the example of the denoising 
algorithm JTJ , we recall the concept of the BM3D modeling. 
The denoising algorithm can be split into three steps. 

1) Analysis. Similar image blocks are collected in groups. 
Blocks in each group are stacked together to form 3-D 
data arrays, which are decorrelated using an invertible 
3D transform. 

2) Processing. The obtained 3-D group spectra are filtered 
by hard thresholding. 

3) Synthesis. The filtered spectra are inverted, providing 
estimates for each block in the group. These blockwise 
estimates are returned to their original positions and the 
final image reconstruction is calculated as a weighted 
average of all the obtained blockwise estimates. 

The blocking imposes a localization of the image on small 
pieces where simpler models may fit the observations. It 
has been demonstrated that a higher sparsity of the signal 
representation and a lower complexity of the model can be 
achieved using joint 3D groupwise instead of 2D blockwise 
transforms. This joint 3D transform dramatically improves the 
effectiveness of image spectrum approximation. 

The total number of groupwise spectrum elements is much 
larger than the image size, and we arrive to an overcomplete or 
redundant data approximation. This redundancy is important 
for effectiveness of the BM3D modeling. 

Our target is to give a strict frame interpretation of the 
analysis and synthesis operations in BM3D. 

A. Matrix representation of analysis and synthesis operations 

Let Y be a y/N x y/N square matrix representing the image 
data and y be the corresponding 1* -vector built from the 
columns of Y. To each y/Nu X y/Nu square image block 
we assign unique index equal to the index of its upper-left 
corner element (pixel) in y. We denote a vector of elements 
of j-th block Yj by yj and define Pj as an Nu x N matrix 
of indicators [0,1] showing which elements of y belong to 
the j-th block, so that yj = Pjy. For the sake of a notation 
simplicity, we assume that the number of blocks in each group 
is fixed and equal to K. Let J r = {j r .i, jV.A'} be the set 
of indices of the blocks in the r-th group, then grouping is 
completely defined by the set J = {J r : r = 1, ...,R}, where 
R is a total number of the groups. It is assumed that for 
each pixel there is at least one block containing the pixel and 
entering in some group. 



The particular form of the 3-D decorrelating transform 
constitutes an important part of the BM3D modeling. It is 
constructed as a separable combination of 2-D intrablock 
and 1-D interblock transforms. The 2-D transform, in turn, 
is typically implemented as a separable combination of 1-D 
transforms. Let D 2 and Di be \/Nbi x \fNbi and K x K 
size matrices representing respectively 1-D interblock and 1- 
D intrablock transforms. Then the separable 2-D transform for 
the block Yj is given by the formula 

0, = D 2 Y,D^. 

The vectorization of this formula using the Kronecker matrix 
product ® gives 

e j = (D 2( g)D 2 ) - yj , 

where 6j , yj 6 R* 61 are the vectors corresponding to the 
matrices ©j and Yj, respectively. To obtain the 3-D spectrum 
of the r-th group we form the Nu x K matrix of the vectorized 
spectrums [6j r 1 , 6j r 2 , 9j r K ] and apply the 1-D interblock 
transform to each row of this matrix 



]r,2 ! •••! "3r,K 



Performing vectorization again, we express the 3-D group 
spectrum coefficients in a compact form: 



E 



ij <g> [(D 2 <x>D 2 ) ■ yj] 



(E, eJr d ^[( D 2® D 2) p .]) -y. 



where ui r is the columnwise vectorized matrix fi r and dj is 
the j-th column of Di. Finally, denoting 



E 



jeJ, 



dj ® [(D 2 ®D 2 ) P 



(3) 



we express the joint 3D groupwise spectrum uj 
[u>l,...,u> T R ] T 
form 



l M of the image Y in the vector-matrix 



'it 



(4) 



The matrix <fr defined by the formulas ©-dUl gives an explicit 
representation of the BM3D analysis operation. 

The synthesis matrix is derived similarly. First, the inverse 
3-D transform is applied to each group spectrum u> r and then 
obtained block estimates are returned to their original positions 
by Pj, j £ J r . The estimate obtained from the r-th group 
spectrum is expressed as ^ r uj r , where 



E d J' 



Py(D 2 ®D 2 ) J 



(5) 



jeJ r 

is an N x Nm matrix. 

The final image estimate is defined as the weighted mean 
of the groupwise estimates using weights g r > 0. Hence the 
synthesis operation has the form 



where 



*u; = W 1 • [g x g R * R ] • w, (6) 
W =E^EpJP J (7) 
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normalizes the weighted mean. W is a diagonal matrix, since 
all products PjPj are diagonal matrices. The m-th diagonal 
element of P JPj is 1 if the m-th pixel of y belongs to the j-th 
block, otherwise it is 0. Thus, the m-th diagonal elements of 
the matrix-sum J2jei P J P i indicates the number of blocks 
in the r-th group containing m-th pixel. 

The matrix VE* defined by the formulas (|5]l-(|2]l gives the 
matrix representation of the BM3D synthesis operation. 

B. Frame interpretation 

Proposition 1: The following equations hold for the matri- 
ces $ and defined by (@ and (0: 

* T • * = E E p J p i > °> w 

r jelr 

* ■ * T = 9r E PjPiW- 2 > 0, (9) 

r jelr 

* -<f> = I NxN . (10) 

The proof is presented in Appendix [A] 

It follows from Proposition Q] that rows of $ constitute a 
frame {</>„} in R . Indeed, let us verify the frame inequality. 
Using the analysis formula u> = <I?y we obtain 

J2\(K,y)\ 2 = " T " = 

n 

If a and b are respectively minimum and maximum values of 

the diagonal matrix ^2 r J2jei P J P J> ^ en ^ or an y y ^ ^ N 
holds the frame inequality 

a -iiyii 2 <£i(0„,y}i 2 <Miyii 2 - < 12 > 

n 

The frame {</>„} is not tight because a ^ b. This follows 
from the fact that the elements on the diagonal of matrix 
X^rSjei PjPj cou nt the number of blocks containing a 
given pixel. These values are different for different pixels, 
since pixels from the blocks possessing higher similarity to 
other blocks participate in a larger number of groups. 

Similarly, using (O we can show that columns of 
constitute a non-tight frame {■»/>„}. From equation (TTOt it 
follows that {</>„} is dual to {ip n }- In general {<fi n } is an 
alternative dual and becomes canonical dual only when all 
weights g r are equal. 

We would like to emphasize that since groups and weights 
are selected data adaptively, the constructed frames are also 
data adaptive. 

The presented frame interpretation allows to extend the 
scope of the BM3D modeling to the modern variational image 
reconstruction techniques. 

III. Variational image deblurring 

The frame based variational image reconstruction problem 
allows two different formulations depending on what kind 
of image modeling, analysis or synthesis is used J4). In 
the analysis formulation the relation between the image and 



spectrum variables is given by the analysis equation u> = $y. 
The problem is formalized as a constrained optimization: 
1 2 

(a>,y) = argmin{^-2 ||z - Ay|| 2 + r • ||w|| p \u = *y}, 

(13) 

where ||-|| is the standard notation of the Z p -norm. 

In the synthesis formulation the relation is given by the 
synthesis equation y = S&u;, leading to the constrained opti- 
mization: 

1 2 

(u>,y) = argminj— -j z - Ay|| 2 +t • w |y = 

(14) 

These problems have equivalent unconstrained forms in which 
they usually encounter in literature. To obtain them it is enough 
to eliminate u> and y respectively from ( fT3l and ( TBI . The 
analysis problem is then formulated as the minimization in 
the image domain 

1 2 

y = argminj^-llz - Ay|| 2 +r ■ ||*y|| p }. (15) 

Similarly, the synthesis problem is formulated as the mini- 
mization in the spectrum domain 

1 2 

Q = argminl-^Hz - A*w|| 2 + r • u> }. (16) 

Despite of the algebraic similarity, the analysis and synthesis 
formulations generally lead to different solutions. A detailed 
discussion of the nontrivial connections between the analysis 
and synthesis formulations can be found in Q. 

The problems (fT3T >-dT6b and the corresponding solution 
techniques recently become a subject of an intensive study. 
In particular, several algorithms have been suggested for 
the convex Zi-norm penalty. These algorithms sharing many 
common ideas are known under different names such as split 
Bregman iterations [8|, iterative shrinkage algorithms ||9l , 
alternating direction method of multipliers [10], majorization- 
minimization algorithms ifTTl . In this paper similar to lfl"2l 
we confine ourself to the Augmented Langrangian (AL) tech- 
nique, using it as a simple and efficient tool for an explicit 
derivation of the reconstruction algorithms. This AL technique, 
introduced independently by Hestenes [13| and Powell lfl4l is 
now widely used for minimization of convex functionals under 
linear equality constraints. 

A. Analysis-based reconstruction 

The AL criterion for the analysis formulation (TPTt takes the 
form: 

L a (y,w,A) = ^ ||z - Ay|£ + t ■ |M| p + 

J-|| w -Sy|£ + i<w-#y,A>X17) 

2 1 7 

where A is a vector of the Lagrange multipliers, 7 > is a 
parameter and the subscript 'a' indicates the analysis formu- 
lation. The saddle problem associated with the Lagrangian L a 
provides the solution of the constrained optimization problem 

<H3. 

Finding the saddle point requires minimization of L a with 
respect to the variables y, ui and maximization with respect to 
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A. A common practical approach is to find the saddle point by 
performing alternating optimization. Applied to ( flTt it results 
in the following iterative scheme: 
Repeat for t = 0, 1, ... 



y t+ i = argminL a (y,oj t , A t ) , (18) 
y 

w t +i = argmini a (y t+ i,u;,A t ) , (19) 

A t+1 = A 4 + ■ (o> t+ i - *y t+1 ) , (20) 

until convergence. 
Here maximization with respect to A is produced as a step 
(|20T > in the direction of the gradient Vx£ a , with a step-size 
/3 > 0. The convergence of the scheme (fT~8b-(l20l> is studied in 

E). 

Minimization with respect to y. Since L a is quadratic with 
respect to y the optimal solution is defined by the linear 
equation 

^A T A+i* T *J .y =-LA T z + i* T (u; + A). (21) 



7 
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We denote by Y a (u, A) the operator giving the solution of 

<ED. 

Minimization with respect to u>. Regrouping the terms in L a 
we arrive to the following formula 



L a (y,o;,A) = — ||z-Ay| 



2a 2 

±\\ U -(*y-\)\\l 

2 7 



H\ P + 
- — Ml 



Since the first and the last terms do not depend on ui, the 
problem is reduced to the optimization 

1 2 

G) = argmiriT • a? + — w - (#y - A) 2 . (22) 

w p 27 

For p < 1, the ^-norm is non-differentiable which makes 
optimization on u) non-trivial. Nevertheless, for p = and 
P = 1 there are well known analytical solutions. 

Let us denote b = <I>y — A, then d22t takes the form 



oj = argnimr • ||w 



i||u;-b||*,w,bel 



(23) 



Depending on the used norm the solution of (123t is given either 
by the hard or soft thresholding according to the formula: 



u> = Th T (b) = 

J 1h* o/t (b) = sign (b) o max (|b| - r, 0) , p=l 
\ £h^(b) = b°l(|b|>v/27), p = 0. 



(24) 



Here all vector operations are elementwise, and 'o' stands for 
the elementwise product of two vectors. We use Th T (b) as 
a generic notation for the thresholding operator. Note, that 
for a given r the thresholding levels for the hard and soft 
thresholdings are calculated differently. 

Applying the general formula (l24l to ( 1221 we obtain the 
solution in the form 

Q = Th T7 (*y - A) . (25) 

Following (fT8l-(|20]i and using d2TT i and d25l l we define 
the analysis-based iterative algorithm which is presented in 



input: z,A,y init 
initialization: 

using yi n i t construct operators $ and 3> T 
set: y ,u;o, A 
t = 
repeat 

A t+ i = A t + /3 • (ut+i - *y t +i) 
t = t + l 

until convergence. 

Fig. 1. Analysis-based deblurring algorithm 



Figure [T] In each iteration it first updates the image estimate 
using the linear filtering ( f2Tb . Then, the difference between 
the spectrum €>y f and A t is thresholded, what corresponds 
to the optimization with respect to u). Finally, the Lagrange 
multipliers are updated in the direction of the gradient u>t+i — 
<fry (+1 . Process is iterated until some convergence criteria is 
satisfied. Particularly, the iterations can be stopped as soon as 
the difference between consecutive estimates becomes small 
enough. 

B. Synthesis-based reconstruction 

The AL criterion for the synthesis formulation (fT4l) takes 
form: 



L 8 (y,o>,A) = 



1 

2^2 



|z - Ay I 



*u;|| 



„., _..„. , -<y-*u>,A). (26) 

2 7 7 

In L s , as opposed to L a , the spectrum variable u> enters the 
quadratic term with a matrix factor \P. It makes the thresh- 
olding formula d24b inapplicable for minimizing i s (y, u>, A) 
with respect to u>. One option is to apply one of the iterative 
shrinkage methods O, but we prefer to follow a different 
approach which leads to a simpler solution. We modify d26l i by 
introducing a splitting variable u 6 R M , used as an auxiliary 
estimate of the spectrum lu. The modified AL takes the form: 

1 2 

Lg(y,w,A,u) = — ||z-Ay||2+r-||u;|| p + 

27 ny 112 
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1 



oj - u 



The corresponding saddle point problem is 

arg min max L s (y, uj, A, u) , 

y,u,u A 



*u, A) + 

(27) 

(28) 



where optimization with respect to the splitting variable u is 
required. 

1 2 

With a small enough £ > penalization by — ||u> — u|| 2 

,2 " ^ 



results in \\uj — ull 



what makes the problem 



equivalent to the saddle problem for 
case we seek for the solution of 



As in the analysis 
by the alternating 
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input: z,A,y init 
initialization: 

using y; n it construct operators & and 1 3/ T 
set: y ,Wo, A ,u 
t = 
repeat 

y*+i=£ (ui,w t) At) 

u t+ i=(7 s (yt,w t+1) A t ) 

w t+ i=Tf) T ^ (u t+ i) 

A t +i = A t + /3 • (y t+ i-*u t+ i) 

t = t + l 
until convergence. 

Fig. 2. Synthesis-based deblurring algorithm 



optimization of £ s (y, u>, A, u) with respect to the variables 
y, u>, u and A. 

Minimization with respect to y is given by the solution of 
the linear equation 



l NxN 



y = — A z H — (*u - A) . (29) 



Minimization with respect to u satisfies the linear equation 



-* T * + Jimkm) • u = -V T (y + A) + -U. (30) 

7^/7 £ 

Minimization with respect to u, thanks to the splitting 
variable u, can be obtained by the thresholding d24b with the 
parameter t£: 

u> = Th T? (u). (31) 

We denote by Y s (u, u>, A) and C/ S (y,w, A) the operators giving 
the solutions of (|29j and (1301 

Using (1291 -OTTl we define the synthesis-based iterative 
deblurring algorithm which is presented in Figure [2] At the 
first two steps the estimates for the image y t and the splitting 
variable u f are updated by solving d29l and (f30b . Then, the 
splitting variable u t+ i is thresholded reducing the complexity 
of the spectrum estimate u>. Finally, the Lagrange multipliers 
are updated in the direction of the gradient yt+i— \Put+i. 
Process is iterated until some convergence criteria is satisfied. 

IV. Decoupling of blur inversion and denoising 

Above we considered algorithms based on the minimization 
of a single objective function. In this section we present an 
alternative approach based on formulation of the deblurring 
as a Nash equlibrium problem for two objective functions. 
This approach allows to split the deblurring problem into two 
subproblems: a blur inversion and denoising, which are then 
solved sequentially. Such a decoupling has several advantages: 

1) The decoupled algorithms are simpler in design and 
parameter selection; 

2) The blur inversion can be implemented efficiently using 
Fast Fourier Transform (FFT); 

3) Various denoising algorithms can be used in this scheme 
selected independently with respect to deblurring; 



4) In many cases decoupled algorithms demonstrate better 
performance than the algorithms where deblurring and 
denoising are performed jointly. 

Examples of the decoupled deblurring can be found in 
works (21, lfl5l . lfl6l and ifTTl . where the regularized inverse 
is followed by different types of filtering (wavelet, shape- 
adaptive DCT, BM3D, pyramidal). An interesting development 
of this technique is demonstrated in [18| where an iterative 
algorithm is derived by alternating optimization of multiple 
objective functions. 

A. Deblurring as a Nash equilibrium problem 

Let us formulate the deblurring problem as the following 
constrained optimization: 



z — Ay||j subjectto ||y — \&u>*|| 2 <£i, 



arg mm 

: argminr • subjectto ||cj — * 1 1 2 



<e 2 , 



(32) 

where E\ , £2 > 0. This problem can be replaced by the 
equivalent unconstrained one: 



y* = argminL inv (y,W*) 
y 

uj* = argmini den (y*,w) 



where 



r(y,«) - ^ 



£den(y, W) 



Ayll^lly- 



(33) 

2,(34) 
(35) 



and 7, £ are constants selected correspondingly to the values 

of e 1 ,e 2 - 

In terms of the game theory the problem d33l can be inter- 
preted as a game of two players identified, respectively, with 
two variables y and us fl9l .l20l. An interaction between the 
players is noncooperative because minimization of L mv (y,uj) 
with respect to y in general results in increase of Lden(y,^) 
and minimization of £den(y,^) with respect to u> increases 
Lj nv (y, <jj). The equilibrium of this game called Nash equilib- 
rium defines the fixed point (y*,o>*) of the optimization. For 
p = 1, problem d33l is convex. 

The objective functions L lm and L^ n allow the following in- 
terpretation. In L mv the fidelity term — ^ ||z — Ay^ evaluates 
the divergency between the observation z and its prediction 
Ay. This fidelity is penalized by the norm ||y — ^w|L 
defining a difference between y and its prediction *&u> through 

uj. The term — - ||a> — 3>y||2 in ^den evaluates a difference 

between the spectrum u> and the spectrum prediction $y 
obtained from y. The error between u> and €>y is penalized 
by the norm ||u>|| p . 

Hence the Nash equilibrium provides a balance between 
the fit of the reconstruction y to the observation z and the 
complexity of the model ||w|L. This can be contrasted with 
the analysis and synthesis-based problem formulations where 
the balance is provided within a single criterion. As we 
demonstrate later the form of the balance plays an essential 
role in the reconstructions with non-tight frames. 
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input: z,A,y init 
initialization: 

using y; n it construct operators 3? and ^ 
set: y ,w = *y 
t = 
repeat 

Deblurring: 



y*+i 



r A T A + -I 



^A T z + ±*u; t 



Deno/s/ng: 
^t+l=1f)r 5 (*y«+i) 

t = t + l 

until convergence. 

Fig. 3. IDD-BM3D - Iterative Decoupled Deblurring BM3D algorithm 



B. IDD-BM3D algorithm 

To solve ( f33l > we consider the following iterative procedure: 
f yt+i = argminL inv (y,u; t ) 



y . T , .,4 = 0,1,... 

= axgmmLden (yt+i,^) 



(36) 



The iterative algorithm d36*i > models the selfish behavior, where 
each variable minimizes only its own objective function. These 
iterations converge to the fixed point (y*,u>*) of d33l l, the 
corresponding result is formulated in Section [V] 

Minimization of L lm with respect to y is given by the 
solution of the linear equation 



— A 1 A + -I ■ y =^A J z + 
a A 1 ) a 7 



(37) 



This step performs regularized inversion of the blur operator. 

The minimization of Lden with respect to u> is obtained by 
thresholding with the threshold parameter r£: 



u = Th T5 (*y) . 



(38) 



Thus, in (|36l l the blur inversion and the denoising steps are 
fully decoupled. 

The algorithm based on (l36*l l is presented in Figure [3] 
We call this algorithm Iterative Decoupled Deblurring BM3D 
(IDD-BM3D)0 

V. Convergence 

A. Analysis and synthesis-based algorithms 

The main motivation of the AL technique is to replace a 
constrained optimization with a simpler saddle-point problem. 
The equivalence of these two problems is not a given fact. 
The classical results stating equivalence are formulated for the 
convex and differentiable functions ||2T1| . Since Z p -norms with 
p < 1 are non-differentiable these results are inapplicable. 
Nevertheless, for the Zi-norm the equivalence can be shown, 

'We wish to note that IDD-BM3D is similar but not identical to our 
Augmented Lagrangian BM3D deblurring (AL-BM3D-DEB) algorithm pre- 
sented earlier in |5|. The AL-BM3D-DEB algorithm is derived from the 
analysis-based formulation (T7J. The regularized inverse step 12 1 1 in AL- 
BM3D-DEB is replaced by the inverse (29} obtained from the synthesis- 
based formulation j26\ . In [5 ] this replacement is treated as an approximation 
and is not mathematically rigorous. The presence of the Lagrange multipliers 
discriminates the AL-BM3D-DEB algorithm from the IDD-BM3D. 



provided that the constraints in the problem are linear. In 
the recent paper l22ll the equivalence statement is proved for 
the total variation penalty. This proof remains valid for any 
convex and non-differentiable penalties, in particularly for the 
Zi-norm based penalties. The equivalence result is formulated 
as following: 

(y, a>) is a solution of the analysis or synthesis problems 
if and only if there exist a saddle-point of the corresponding 
ALs. 

Practically it means that the saddle-point of the AL opti- 
mization can be used in order to obtain the solutions of the 
considered optimization problems. 

The convergence properties for the analysis and synthesis- 
based algorithms are formulated in the following proposition. 

Proposition 2: 

(a) If there exists a saddle point (y*, ui* , A*) of L a (y, u;, A) 
( 1771 ), then y t — !• y*.u; t ^ u}*,\-^ A*. 

(b) If there exists a saddle point 
(y*,u;*,u*,A*) of L s (y,u>,u,A) (ETJ), then 
y«-> y*,w t -> ^*,u t ^ u*, A t -> A*. 

On the other hand, if no such saddle point exists, then at 
least one of the sequences {y t } or {A t } must be unbounded. 
The proof is given in Appendix IB"1 

B. IDD-BM3D algorithm 

Proposition 3: For any set of parameters er, r, 7, £ the 
sequence (yt, w t ) generated by the IDD-BM3D algorithm with 
equal group weights g r , converges to the fixed point (y*,w*) 
defined by the equations ( 1331 ), the fixed point exists. 

The proof of the proposition is given in Appendix [B] It is 
not required that the fixed point is unique. Depending on a 
starting point (yo, Wo) the limit point of the algorithm can be 
different but should satisfy the fixed point equations. 

VI. Implementation 

Grouping and frame operators. To build the groups, we use 
the block-matching procedure from HI and apply it to the 
image reconstructed by the BM3DDEB deblurring algorithm 
0. The found locations of the similar blocks constitute the 
set J that is necessary to construct the analysis and synthesis 
frames. Multiplications against the matrices $,# T ,* and 
\f T are calculated efficiently since all of them involve only 
groupwise separable 3-D transformations of the data (possibly 
with some averaging of the estimates). In our experiments the 
3-D transform is performed by first applying the 2-D discrete 
sine transform (DST) to each block in the group followed by 
the 1-D Haar transform applied along the third dimension of 
the group. The image block size is 4 x 4, and the number of 
blocks in the group is 8. 

Choice of the group weights. Since image blocks are 
overlapping, for each pixel we obtain several estimates. The 
weighted averaging can be used to improve the final aggre- 
gated estimate. For the one-step (non-iterative) algorithms the 
weights can be adaptively selected so to minimize the variance 
of the final aggregated estimate, based on the variance of 
each of the estimates (e.g. l23l . 0], 0). In the considered 
iterative algorithms the influence of the weights on the final 
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estimate is complex, and deriving a formula for the optimal 
weights is rather involved. Instead, following the idea of 
the sparse representations, we suggest giving the preference 
to the estimates obtained from the sparser groups. In our 
implementations we use weights inversely proportional to the 
number of significant spectrum coefficients of the groups 
g r = 1/ ||Th e (av)|| , where significant coefficients are found 
by the hard thresholding of the group spectra using a small 
threshold e. 

The grouping and the adaptive group weights are calculated 
only once, using the initial image estimate yi n it and remain 
unchanged through the subsequent iterations. 

Choice of the regularization parameters. The parameters 
t, 7, £ are optimized to provide best reconstruction quality. 
Optimization has been performed separately for each algo- 
rithm and each deblurring scenario. The parameter f3 is always 
set to 1. 

Initialization. We experimentally confirmed the convergence 
to an asymptotic solution that is independent of the initial- 
ization yo and wo- Nevertheless, initialization with a better 
estimate, for example with the reconstruction obtained by 
BM3DDEB (which we also use to define grouping) results 
in a much faster convergence. 

Solution of the large-scale linear equations. All proposed 
algorithms contain steps involving solution of large-scale 
linear equations. For a circular shift-invariant blur operator, the 
solution of the equations ((29) and rt37l can be calculated in the 
Fourier domain using the FFT The more complex equations 
(|2TT > and ( f3Qb are solved using the conjugate gradient method. 
The conjugate gradient method allows avoiding explicit calcu- 
lations of the matrices $ T $ and ^ T, ii, since it requires only 
evaluating products of these matrices against vectors. 

Practical considerations. The two steps of the IDD-BM3D 
algorithm can be merged into a single one 



yt+i = F~ 



'T* (h)o^(z) + ^(*Sf) r€ (*y t ))" 



+ — 
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where the analysis-thresholding-synthesis operation 
fy1t) T £ (3?yt) can be calculated groupwise without need 
to obtain the whole spectrum id t explicitly. Here h denotes 
the vectorized blurring kernel corresponding to the blur 
operator A, and 'o' stands for the elementwise product of 
two vectors. The operator T (■) reshapes the input vector into 
a 2-D array, performs 2-D FFT and vectorizes the obtained 
result. J 7-1 (•) works analogously, performing inverse FFT. 

Complexity. Application of the frame operators is the most 
computationally expensive part of the proposed algorithms. 
However, due to their specific structure, the complexity of the 
frame operators 3? and \& is growing only linearly with respect 
to the number of the pixels in the image. To give an estimate 
of the complexity of the IDD-BM3D algorithm, we mention 
that, on a 256 x 256 image, one iteration takes about 0.35 
seconds, and about 50 iterations are typically sufficient. This 
timing has been done on dual core 2.6 GHz processor for 
an implementation where the computationally most intensive 
parts have been written in C++. 



Scenario 


PSF 


(7 


1 


1/ (1 + x'f + xjj) , xi,X2 = —7, ■ ■■ ,7 


2 


2 


1/ (1 + xi + x'fj , X!,X2 = -7, . . . , 7 


8 


3 


9x9 uniform 


ks 0.3 


4 


[1 4 6 4 l] 1 ' [1 4 6 4 1] /256 


49 


5 


Gaussian with std =1.6 


4 


6 


Gaussian with std = 0.4 


64 



TABLE I 

Blur PSF and noise variance used in each scenario. 




on 
z 
<n 



hard thresholding, initialized with BM3DDEB 
hard thresholding, initialized with 
soft thresholding, initialized with BM3DDEB 
soft thresholding, initialized with 



50 



100 
Iterations 



150 



200 



Fig. 4. Change of the ISNR with iterations for the different setups of the 
IDD-BM3D algorithm. Deblurring of Cameraman image, scenario 3. 



VII. Experiments 

We consider six deblurring scenarios used as the bench- 
marks in many publications (e.g., lfT7l and flU). The blur point 
spread function (PSF) h (xi,X2) and the variance of the noise 
cr 2 for each scenario are summarized in Table U PSFs are 
normalized so that ^ h = 1. Each of the scenarios was tested 
with the four standard images: Cameraman, Lena, House and 
Barbara. 



A. Experiment 1 - comparison of the proposed algorithms 

All three proposed algorithms, namely: analysis-based, 
synthesis-based and IDD-BM3D are evaluated in the scheme 
with the soft thresholding and unit group weights (g r — 1). 
Additionally, the IDD-BM3D algorithm is tested with the 
adaptive group weights (g r — 1/ ||Xh e (u; r )|| ) using the soft 
and hard thresholdings. 

In Table [TT] we present improvement of signal-to-noise ratio 
(ISNR) values achieved by each algorithm for the Cameraman 
image. From these values we can conclude that the synthesis- 
based algorithm performs essentially worse than the IDD- 
BM3D algorithm, with the analysis-based algorithm being in- 
between. We can also see that the adaptive weights indeed 
provide a noticeable restoration improvement. Finally, com- 
paring the last two rows, we conclude that hard thresholding 
enables better results than the soft thresholding, and combined 
with the adaptive weights it provides the best results among 
the considered algorithms. 

Convergence properties of the IDD-BM3D algorithm are 
demonstrated in Figure |4] 



The experiments with the IDD-BM3D algorithm can be 
reproduced using the Matlab program available as a part of 
the BM3D package^. 

B. Experiment 2 - comparison with the state of the art 

Table [III] presents a comparison of the IDD-BM3D algo- 
rithm versus a number of algorithms including the current state 
of the art. The ISNR values for ForWaRD |24), SV-GSM Q21, 
SA-DCT HH and BM3DDEB [2 1 are taken from our previous 
paper 0, while the results for LO-AbS [25 1, TVMM ifTTl . 
CGMK ll26ll are obtained by the software available online. 
We use the default parameters suggested by the authors of 
the algorithms. The IDD-BM3D algorithm in this comparison 
employs the hard thresholding and the adaptive weights. 

The proposed IDD-BM3D algorithm provides the best re- 
sults with significant advantage over closest competitors. Par- 
ticularly interesting is the comparison against the BM3DDEB 
algorithm. BM3DDEB is a two-stage non-iterative algorithm. 
On the first stage it utilizes the BM3D image modeling to 
obtain the initial estimate, which is then used on the second 
stage for an empirical Wiener filtering. Better performance 
of the IDD-BM3D algorithm demonstrates that considered 
decoupled formulation ( f33l > enables more effective exploit- 
ing of the BM3D-modeling than the two-stage approach of 
BM3DDEB. 

The visual quality of some of the restored images can be 
evaluated from Figures [5] and [6] where for a comparison we 
show results by the closest competitors |26|, |25| and |2|. One 
can see that the proposed algorithm is able to suppress the 
ringing artifacts better than BM3DDEB and provides sharper 
image edges. This latter effect is achieved in particular due 
to the smaller block size used in IDD-BM3D compared to 
BM3DDEB. 

VIII. Discussion 

In the experiments of the previous section we observed a 
clear advantage of the IDD-BM3D algorithm over the analysis- 
based one. This result is rather surprising, since in the case 
of the tight frames the IDD-BM3D and the analysis-based 
algorithms are almost identical. 

Indeed, if we assume that {</>„} is a tight frame and require 
that all group weights will be equal, then = al and 

\I> = (<& T <I?) <& T = a _1 $ T . Substituting these expressions 
into equation (f2TT > of the analysis-based algorithm we obtain 

( i,A T A + -iV y = Aa t z + (« + A) . 

Comparing it with the equation d37b we see that up to the pres- 
ence of the Lagrange multipliers the analysis-based algorithm 
is identical to the IDD-BM3D algorithm. This observation 
rises a question: what makes the algorithms behave differently 
when the frame is not tight? 

To find an answer, let us look again at the equation ( BTT i. Its 
solution requires inversion of the matrix ^jA T A + i$ T 4>, 
whoes condition number depends not only on the properties 

1http://www.cs.tot.firfoi/GCF- BM3D 



of the blur operator but also on the properties of the frame. 
In the case of the non-tight analysis BM3D-frame, 3> T <I? is a 
diagonal matrix, its entries are defined by the data grouping 
and count number of times each pixel appears in different 
groups. Experiments demonstrate that the variation of these 
entries can be very large (up to hundreds times). The large 
differences in magnitude of the diagonal elements of €> T $ 
make the matrix -lA T A+-$ r $ ill-conditioned and result in 
degradation of image reconstruction compared to IDD-BM3D. 

Presence of the matrix 3> T <I> in the reconstruction formulas 
is inevitable as long as one uses criterion containing norms 
both for the image and spectrum domain. Formulation based 
on the Nash equilibrium allows to overcome this problem and 
have norms only from one domain in each criterion. 

IX. Conclusions 

The frame based formulation opens new perspectives for the 
use of BM3D modeling within the variational reconstruction 
techniques. The developed deblurring algorithm demonstrates 
state-of-the-art performance, confirming a valuable potential of 
BM3D-frames as an advanced image modeling tool. For non- 
tight frames, we argue the validity of image reconstruction 
by minimizing a single objective function and propose an 
alternative formulation, based on Nash equilibrium problem. 

Appendix A 

A. Proof of Proposition [7] 

The proof is based on use of the following Rronecker matrix 
product formulas. 

If A is an m x n matrix and B is a p x q matrix, then the 
Rronecker product A ® B is the rap x nq block matrix and 

(A®B)(C®D) = AC ® BD, 
(A<g)B) T = A T (g)B T , 
(A^B)- 1 = A -1 (SB -1 . 

Also, matrix equation AXB = C can be vectorized column- 
wise with respect to X and C as following 

(B T ® A)vect (X) = vect (C) . 

To simplify notation we denote G = (Di ® Di), Then the 
formula © from Proposition Q] is proved as following 

* T * = *r *r = 

EE E( d J®pjG T )(d,i^GP,o = 
EE E (djd,o ® (pjg-gp,) = 

r jeJrj'eJ-r 

EE E^/pj-r-p/ EE^tv 

r jeJrj'EJr r j£j r 




Fig. 6. Deblurring of the Lena image, scenario 2. From left to right and from top to bottom are presented zoomed fragments of the following images: 
original, blurred noisy, reconstructed by CGMK (26) (ISNR 5.37), LO-AbS (25] (ISNR 5.71), DEB-BM3D |2| (ISNR 6.53) and by proposed IDD-BM3D 
method (ISNR 6.61). 
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Method 


Scenario 


| Thresh. | Weights g r 


1 | 2 | 3 1 4 | 5 | 6 


Cameraman (256x256) 


BSNR 






31.87 


25.85 


40.00 


18.53 


29.19 


17.76 


Input PSNR 






22.23 


22.16 


20.76 


24.62 


23.36 


29.82 


Synthesis 


soft 


unit 


6.30 


4.60 


7.88 


2.06 


2.98 


2.84 


Analysis 


soft 


unit 


7.88 


5.75 


9.22 


3.00 


3.67 


3.92 


IDD-BM3D 


soft 


unit 


8.17 


6.17 


9.38 


3.17 


3.83 


4.12 


IDD-BM3D 


soft 


adaptive 


8.41 


6.41 


9.59 


3.38 


3.98 


4.14 


IDD-BM3D 


hard 


adaptive 


8.85 


7.12 


10.45 


3.98 


4.31 


4.89 



TABLE II 

Comparison of the output ISNR [dB] of the proposed deblurring algorithms. Row corresponding to "Input PSNR" contain PSNR 
[dB] of the input blurry images). Blurred signal-to-noise ratio (BSNR) is defined as 10log\o (var (Ay) /Na 2 ), where var() is the 

variance. 





[ Scenario 


[ Scenario 




1 | 2 | 3 | 4 1 5 1 6 


1 1 2 1 3 | 4 I 5 1 6 


Method 


| Cameraman (256x256) 


| House (256x256) 


BSNR 


31.87 


25.85 


40.00 


18.53 


29.19 


17.76 


29.16 


23.14 


40.00 


15.99 


26.61 


15.15 


Input PSNR 


22.23 


22.16 


20.76 


24.62 


23.36 


29.82 


25.61 


25.46 


24.11 


28.06 


27.81 


29.98 


ForWaRD [24] 


6.76 


5.08 


7.34 


2.40 


3.14 


3.92 


7.35 


6.03 


9.56 


3.19 


3.85 


5.52 


SV-GSM [17J 


7.45 


5.55 


7.33 


2.73 


3.25 


4.19 


8.64 


7.03 


9.04 


4.30 


4.11 


6.02 


SA-DCT 1 16 1 


8.11 


6.33 


8.55 


3.37 


3.72 


4.71 


9.02 


7.7 '4 


10.50 


4.99 


4.65 


5.96 


BM3DDEB |2| 


8.19 


6.40 


8.34 


3.34 


3.73 


4.70 


9.32 


8.14 


10.85 


5.13 


4.56 


7.21 


LO-AbS 1 25 1 


7.70 


5.55 


9.10 


2.93 


3.49 


1.77 


8.40 


7.12 


11.06 


4.55 


4.80 


2.15 


TVMM 1 1 1 1 


7.41 


5.17 


8.54 


2.57 


3.36 


1.30 


7.98 


6.57 


10.39 


4.12 


4.54 


2.44 


CGMK [26 ] 


7.80 


5.49 


9.15 


2.80 


3.54 


3.33 


8.31 


6.97 


10.75 


4.48 


4.97 


4.59 


IDD-BM3D 


8.85 


7.12 


10.45 


3.98 


4.31 


4.89 


9.95 


8.55 


12.89 


5.79 


5.74 


7.13 




| Lena (512x512) 


| Barbara (512x512) 


BSNR 


29.89 


23.87 


40.00 


16.47 


27.18 


15.52 


30.81 


24.79 


40.00 


17.35 


28.07 


16.59 


Input PSNR 


27.25 


27.04 


25.84 


28.81 


29.16 


30.03 


23.34 


23.25 


22.49 


24.22 


23.77 


29.78 


ForWaRD |24| 


6.05 


4.90 


6.97 


2.93 


3.50 


5.42 


3.69 


1.87 


4.02 


0.94 


0.98 


3.15 


SV-GSM 1171 














6.85 


3.80 


5.07 


1.94 


1.36 


5.27 


SA-DCT |16| 


7.55 


6.10 


7.79 


4.49 


4.08 


5.84 


5.45 


2.54 


4.79 


1.31 


1.02 


3.83 


BM3DDEB |2| 


7.95 


6.53 


7.97 


4.81 


4.37 


6.40 


7.80 


3.94 


5.86 


1.90 


1.28 


5.80 


LO-AbS [25 




6.66 


5.71 


7.79 


4.09 


4.22 


1.93 


3.51 


1.53 


3.98 


0.73 


0.81 


1.17 


TVMM 1 1 1 




6.36 


4.98 


7.47 


3.52 


3.61 


2.79 


3.10 


1.33 


3.49 


0.41 


0.75 


0.59 


CGMK |26 | 


6.76 


5.37 


7.86 


3.49 


3.93 


4.46 


2.45 


1.34 


3.55 


0.44 


0.81 


0.38 


IDD-BM3D 


7.97 


6.61 


8.91 


4.97 


4.85 


6.34 


7.64 


3.96 


6.05 


1.88 


1.16 


5.45 



TABLE III 

Comparison of the output ISNR [dB] of deconvolution methods (row corresponding to "Input PSNR" contain PSNR [dB] of the 

input blurry images). 



Proof of the formula ((9]): 

(W- 1 -[.gi*i,..., 5 K*fl]) T = 
W- X [si • [si*i,...,5fl* H ] T W- 1 = 

w E";E E<w ( p J p ;) w 1 

r jeJ r j'eJ r 



w E'^E p J p i- 

r j£j r 

The last identity holds since J2 r 9rJ2jeJ PjPj an< ^ ^ 
are diagonal matrices. 



The formula dTOb in Proposition Q] is valid since 





= (W- 1 • [(ft*! 




w 


r 




w 




SPJGT )l 




r \ jeJr 




w 


'E^E E 


( d J^PjG T 




r jeJrj'GJr 




w 


'E^E E 


(djd s ) « (P 




r jeJrj'eJr 




w 


X E^EE 


% (PfP, ) = 




r jeJrfeJr 




w 


^EpJp 


j = Inxn- 
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Appendix B 

A. Proof of Proposition [2] 

Let us consider constrained optimization problem given in 
the following general form 

9 



(a) Comparing the AL with (O we note that / (u) = 
— -77 llz — Ay II? and the equality Cv+Du = b takes the form 
uj — <I>y = 0, where u) corresponds to v and u corresponds 
to y. Thus, C = Ij,/xm an d D = — 

We have two conditions of the theorem to be tested: C 

2 • 



min{/ (u) + V 9j ( Vi ) |Cv + Du = b}, (39) has ful1 column rank and f ^ + H Du lla is strictl y con J ex 



In our case, C = I 



MxM 



has full column rank, IIDul 



where u £W n , v = [vj, vj] T , Vj £ R m * , v £ R m , fh = 
^Tij,b £ R S ,C is of the size (s x m),D is of the size 
(s x m) and /(u) is convex. The AL corresponding to this 
problem is 

9 

L(u,v,A) = /(u) + ^<?i(vj) + 

a ||Cv + Du b||2 + (Cv + Du b, A) (40) 

The link between the main variable u and the auxiliary split- 
ting variable v is given by the linear equation Cv + Du = b. 
If C is the identity matrix, then v = b Du and the conver- 
gence of the corresponding iterative algorithm can be obtained 
from the Eckstein-Bertsekas's theorem ([21], Theorem 8). 
However, if Cv + Du = b is not resolved with respect to 
v then the theorem is not applicable in its original form. 
The techniques exploited in our paper leads to the relations 
between the variables which cannot be resolved with respect 
to v. In order to analyze the convergence of the proposed algo- 
rithm we use a novel formulation of the Eckstein-Bertsekas's 
theorem ||27l adapted to the general linear link between the 
variables v and u. This new Eckstein-Bertsekas's theorem is 
given in the following form (27]. 

Theorem 4: Consider the problem d39l where / and gj are 
closed proper convex functions, C has full column rank and 
/ (u) + HDulla is strictly convex. Let u a £ M m , A € M s be 
arbitrary and /3 > 0. Suppose that there are sequences {er^ } 
and {v t } such that a\ > 0, Vt > and J2 t a\ < oo, J2 t v t < 
oo. Assume that 



v t+ i - argmin |E*=i 9j (vj) + 
+a ||Cv + Dut-bHa + (Cv, A 



u t+ i - argmin 



+ 



V + l 



+a||Cv t+ i +Du-b||2 + 
= \ t + /3 (Cv t+ i + Du t+ i 



;du,a ( 

-b). 



< v t , 



If there exists a saddle point (v*, u*, A*) for L (u, v, A) d40l >, 
then vt — > v*,ut — > u*,X t — > A*. On the other hand, if no 
such a saddle point exists, then at least one of the sequences 
{ut} or {Af} must be unbounded. 

This formulation of the convergence concerns approximate 
solutions on each optimization step, where the parameters a\ 
and Vt controls the accuracy at each step. The finite sums 
J2t a t < °°>Et^ < 00 mean that o\,vt 0, i.e. the 
accuracy should asymptotically improve. 

Armed with this theorem we can proceed to the proof of 
Proposition [2] 



,1 



(* T *u,u). Due to © * T * = W>0, thus ||Du||2 is 

strongly convex and the same holds for — - ||z — Ay)^ + 

|Du|| 2 . Thus, all conditions of the theorem are satisfied and 
the analysis-based algorithm converges to the saddle-point 
of the AL ( TTTb , if it exists. It proves the first part of the 
proposition. 

(b) Comparing the formulation d26l l with (|39l we note that 
/ (u) = — 77 | j sb — Ayllj and the equality Cv + Du = b 

takes the form y — SPu = and u> — u = 0. Assuming v — > 
y\ 

, u —> u) these equations give 



lATxJV 





IMxM 



D 



'NxM 
ImxM 



The matrix C is square triangular with elements of the main 
diagonal equal to 1. It has full column rank. For ||Du|| 2 we 
have ||Du|| 2 — > \\u>\L- Thus ||Du|L is strongly convex and 
the both conditions of the theorem are fulfilled. It proves the 
second part of the proposition. 



B. Proof of Proposition \3\ 

We consider the IDD-BM3D algorithm with soft threshold- 
ing and equal group weights g r = c, c € M + ,r = 1,...,R. 
From ©, ©, (01 and ® follows that $ T * = W and 
* = W _1 * T . 

Each iteration of the IDD-BM3D algorithm consists of two 
steps 

' yt+^M- 1 [£A T z + *uJt] , 
u> t+ i = 1h re (*y t+1 ) , 

where M = ^A T A + I > 0. 



(41) 



(42) 



Introducing the operator O d (u) = #M [^-A T z + tyui 
and denoting q t = &Yt we rewrite d4~TT) in a compact form 

qt+i = O d (w t ) , 
W t +l = ^r? (qt+i) • 

The convergence analysis is based on the technique of 
nonexpansive operators. An operator P : M m — > M m is called 
nonexpansive if for any x, x' £ M. m 

||P(x)-P(x')||2<||x-x'||^. 

It is shown in ll28l (Proposition 3.1) that the soft threshold- 
ing is a nonexpansive operator 



Tf£ o/t (x) - £f)* o/t (x') 
with equality holding only when 



< llx-x'| 



2 ' 



1l) s T ° ft (x) - Tf)^* (x') = x - x'. 
Hence the operator Tf) r ^ (•) in ( 1421 is nonexpansive. 



(43) 
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To prove that the operator O d in (l42l i is also nonexpansive, 
we first notice that 

O d (w) - O d (fa/) = ^M" 1 * (w - a/) . 

To find the norm of the matrix <I>M _1 >t' we evaluate 
its eigenvalues. For the matrix <&M~ 1 \I', the corresponding 
characteristic equation is defined as a determinant of the 
equation 

(<&M _1 W X <& T - AI)v = 0, (44) 

where v £ K M is an eigenvector and A is an eigenvalue. The 
matrix 3>M ^ has the size M x M while its rank is equal 
to N. Thus, M — N eigenvalues of this matrix are equal to 
zero. We wish to show that nonzero eigenvalues of ^M -1 ^ 
coincide with the eigenvalues of the matrix M _1 . 

Let us replace in d44b v by <frv, v £ M. N , and multiply the 
equation ( l44b by W _1 <J> T . Then, this equation takes the form 

W- 1 * T (*M- 1 W _1 * T - AI)*v = 0. (45) 

Multiplication by W _1 $ T in ( f45l > is legitimate because it 
preserves the rank of this system of the linear equations. Since 
W _1 * T * = I, (|45ll takes the form 

(M" 1 - AI)v = 0. (46) 

Here A and v become the eigenvector and eigenvalue for 
the matrix M _1 . The eigenvalues of the matrix M _1 = 
[^■A T A + 1] are positive and take values less than or 
equal to 1, 

The passage from (l44l to (|46| > proves that nonzero eigenval- 
ues of the matrix $M _1 \f are equal to the eigenvalues of the 
matrix M _1 .Thus all eigenvalues of the matrix 3>M ^ are 
nonnegative and take values less than or equal to 1. Hence, 
the matrix norm p (ffrM -1 *) is less than or equal to one, 
and the operator Od is nonexpansive due to the inequality 

||O d (w) - O d (fa/)|| 2 = H^M- 1 * (fa, - fa/)|| 2 

< p (^IVr 1 *) - oj'|| 2 < ll w - w 'll 2 ■ 

Let (y*,u>*) be a fixed point of the equations fiTT i and 
A y*= y* - y*, Afa^ = fa>t - fa>*, Aq t = *Ay. Since Xfj T? 
and Od are nonexpansive operators we have from (l42l i that 
||Aq t+ i|| < ||Afa, t || and ||Ao; f+ i|| < ||Aq t+ i||. It follows 
that ||Afa>t_|_i|| < ||AfaJt|| for Vt. Then, the sequence u>t+i 
lies in a compact region and converging to a limit point, say 
w, liirik^oo \\<*>t k — w *| = ||^ — i- e - a distance from 
this limit point to a fixed point is bounded. By the continuity 
of the operators in d4lT i the same statement holds for the 
sequence y t : at least one limit point exists, denoted as y, 
and a distance between this limit point and a fixed point is 
bounded, Hindoo ||y tfc -y*|| = ||y-y*||. 

Again due to the continuity of the operators in (|4"TT i the 
limit point is a fixed point. Replacing (y*,fa>*) by (y,u>) 
we obtain the convergence of the decoupling algorithm, 
lirrik^oo ||w tfc - u>|| = and Um k ^oo \\yt h ~ 9\\ = 0. It 
proves Proposition |3] 
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